home *** CD-ROM | disk | FTP | other *** search
/ Night Owl 6 / Night Owl's Shareware - PDSI-006 - Night Owl Corp (1990).iso / 039a / d3d.zip / HLPFUN.C < prev    next >
C/C++ Source or Header  |  1991-08-30  |  27KB  |  938 lines
  1. /* HLPFUN: A function for hidden- line                      */
  2. /* elimination, using device-independent pixels.            */
  3. /* This version includes a providion for curved surfaces    */
  4.  
  5. #include <stdio.h>
  6. #include <math.h>
  7. #include <ctype.h>
  8. #include <alloc.h>
  9. #include <conio.h>
  10. #include <string.h>
  11. #include <process.h>
  12. #include "grpack.h"
  13.  
  14. #define  max(x,y)    ((x)>(y)?(x):(y))
  15. #define  min(x,y)    ((x)<(y)?(x):(y))
  16. #define  max3(x,y,z) ((x)>(y)?max(x,z):max(y,z))
  17. #define  min3(x,y,z) ((x)<(y)?min(x,z):min(y,z))
  18. #define  xwhole(x)   ((int)(((x)-xmin)/deltaX))
  19. #define  ywhole(y)   ((int)(((y)-ymin)/deltaY))
  20. #define  xreal(i)    (xmin+((i))*deltaX)
  21.  
  22. #define  M  -1000000.0
  23. #define  NSCREEN  15
  24. #define  BIG   1.e30
  25.  
  26. int abs(int i);
  27.  
  28. static void
  29.    check_kbhit(void),
  30.    add_linesegment(int P, int Q),
  31.    viewing(double x, double y, double z,
  32.             double *pxe, double *pye, double *pze),
  33.    linesegment(double xP, double yP, double zP,
  34.                double xQ, double yQ, double zQ, int k);
  35.  
  36. void
  37.    ermes(char *str),
  38.    coeff(float rho, float theta, float phi),
  39.    init_viewport(void);
  40.  
  41. static int
  42.    counter_clock(int i0, int i1, int i2, int code),
  43.    includesvertex(int h, int i, int j),
  44.     sameside(double xj, double yj, double xl, double yl,
  45.                 double xh, double yh, double xi, double yi),
  46.    allocsize(int divfactor, int elsize);
  47.  
  48. extern int nmax, nface, **pface;
  49.  
  50. extern float xO, yO, zO, xmin, xmax, ymin, ymax;
  51.  
  52. static int ipixmin, ipixmax, ipixleft, ipixright,
  53.    ipix, jpix, jtop, jbot, j_old, jI, topcode[3], *POLY,
  54.    npolyalloc, *vertexconvex, ntrsetalloc, npoly,
  55.    isize=sizeof(int), LOWER[NSCREEN], UPPER[NSCREEN],
  56.    LOW[NSCREEN], UP[NSCREEN], ntrset, maxvertex, jface,
  57.     nTRIANGLE, nTRLIST, iTRIANGLE, iTRLIST, inode, itria,
  58.    *trset;
  59.  
  60. extern double v11, v12, v13, v21, v22, v23, v31, v32, v33, v43,
  61.    PIdiv180;
  62.  
  63. static double d, c1, c2,
  64.     eps = 1e-5, meps = -1e-5, oneplus=1.+1.e-5,
  65.     Xrange, Yrange, Xvp_range, Yvp_range,
  66.    deltaX, deltaY, denom, slope,
  67.    Xleft[3], Xright[3], Yleft[3], Yright[3],
  68.     a, b, c, surfl;
  69.  
  70. struct linsegface {
  71.    int i;
  72.    double a, b, c;
  73.    struct linsegface *next;
  74. } *lsegfacenode(void);
  75.  
  76. struct vertex {
  77.    int inuse;
  78.    float xw, yw, zw;
  79.    double xe, ye, ze;
  80.    struct linsegface *connect;
  81. } *VERTEX, *pvertex;
  82.  
  83. extern struct vertex *p;
  84.  
  85. static struct triangle {
  86.    int A, B, C;
  87.    double a, b, c, h;
  88. } *TRIANGLE, *ptriangle;
  89.  
  90. static struct lseglist {
  91.    double xP, yP, zP, xQ, yQ, zQ;
  92.    int k;
  93.    struct lseglist *next;
  94. } *lsegstart, *auxlseg, *lsegnode(void);
  95.  
  96. struct trlist {
  97.    unsigned int ptria;
  98.    unsigned int next;
  99. } *TRLIST, *pnode;
  100.  
  101. struct {
  102.    int tr_cov;
  103.    float tr_dist;
  104.    unsigned int start;
  105. } SCREEN[NSCREEN][NSCREEN], *pointer;
  106.  
  107. void hlpfun(float rho, float theta, float phi, float surflimit)
  108. {
  109.    int P, Q, ii, vertexnr, i, kk, i0, i1, i2, code, count,
  110.       polygonconvex, loopcount, jtr, j, n, k;
  111.  
  112.    struct linsegface *ptr, *dummy, *dummy0;
  113.  
  114.    double
  115.       Xvp_min=0, Xvp_max, Yvp_min=0, Yvp_max,
  116.       fx, fy, Xcenter, Ycenter, Xvp_center, Yvp_center,
  117.       xP, yP, zP, xQ, yQ, zQ, XP, YP, XQ, YQ,
  118.       Xlft, Xrght, Ylft, Yrght,
  119.       xxP, yyP, zzP, xxQ, yyQ, zzQ;
  120.  
  121.    textXY(40, 40, "Please wait ...");
  122.    Xvp_max = x_max;
  123.    Yvp_max = y_max;
  124.    surfl = surflimit;
  125.    npolyalloc = 200;
  126.    POLY =         (int *)farcalloc((long)npolyalloc, (long)isize);
  127.    vertexconvex = (int *)farcalloc((long)npolyalloc, (long)isize);
  128.    if(POLY == NULL || vertexconvex == NULL) {
  129.       ermes("Memory: polygon");
  130.    }
  131.  
  132.    VERTEX = p;
  133.       /* Now VERTEX points to the same memory area as p does
  134.          in module D3D  */
  135.    
  136.    /* Initialize screen matrix */
  137.    for(ipix=0; ipix<NSCREEN; ipix++)
  138.    for(jpix=0; jpix<NSCREEN; jpix++) {
  139.       pointer = &(SCREEN[ipix][jpix]);
  140.       pointer->tr_cov = -1;
  141.       pointer->tr_dist = BIG;
  142.       pointer->start = -1;
  143.    }
  144.  
  145.    coeff(rho, theta, phi);
  146.  
  147.    maxvertex = nmax - 1;   /*nmax is defined in module D3D  */
  148.    dummy = lsegfacenode();
  149.    for(i=1; i<=maxvertex; i++) 
  150.       VERTEX[i].connect = (VERTEX[i].inuse ? NULL : dummy);
  151.            /* dummy indicates: not in use               */
  152.            /* NULL  indicates: empty list (but in use)  */
  153.  
  154.    /* Compute screen constants */
  155.    Xrange = xmax - xmin;
  156.    Yrange = ymax - ymin;
  157.    Xvp_range = Xvp_max - Xvp_min;
  158.    Yvp_range = Yvp_max - Yvp_min;
  159.    fx = Xvp_range/Xrange;
  160.    fy = Yvp_range/Yrange;
  161.    d = (fx<fy ? fx : fy);
  162.    Xcenter = 0.5*(xmin + xmax);
  163.    Ycenter = 0.5*(ymin + ymax);
  164.    Xvp_center = 0.5*(Xvp_min + Xvp_max);
  165.    Yvp_center = 0.5*(Yvp_min + Yvp_max);
  166.    c1 = Xvp_center - d*Xcenter;
  167.    c2 = Yvp_center - d*Ycenter;
  168.    deltaX = oneplus*Xrange/NSCREEN;
  169.    deltaY = oneplus*Yrange/NSCREEN;
  170.    /* Now we have: Xrange/deltaX < NSCREEN */
  171.  
  172.    iTRLIST = 0;
  173.    iTRIANGLE = 0;
  174.  
  175.    nTRLIST = allocsize(4, sizeof(struct trlist));
  176.     TRLIST  = (struct trlist *)farcalloc((long)nTRLIST,
  177.                         (long)sizeof(struct trlist));
  178.  
  179.    nTRIANGLE = allocsize(3, sizeof(struct triangle));
  180.    TRIANGLE = (struct triangle *)farcalloc((long)nTRIANGLE,
  181.                         (long)sizeof(struct triangle));
  182.  
  183.    if(TRLIST == NULL || TRIANGLE == NULL)
  184.       ermes("Memory: triangles");
  185.  
  186.    for(j=0; j<nface; j++) {
  187.       if(pface[j] == NULL)
  188.          continue;
  189.       n = pface[j][0];
  190.       if(n < 0)
  191.          ermes("Internal error: n < 0");
  192.       if(n == 1)
  193.          ermes("Polygon with only one vertex");
  194.       POLY[0] = pface[j][1];
  195.       npoly = 1;
  196.       for(k=2; k<=n; k++) {
  197.          i = pface[j][k];
  198.          if(npoly == npolyalloc) {
  199.             npolyalloc += (npolyalloc>>2);
  200.             POLY = (int *)farrealloc(POLY,(long)npolyalloc*isize);
  201.             vertexconvex = (int *)
  202.                farrealloc(vertexconvex,(long)npolyalloc*isize);
  203.             if(POLY == NULL || vertexconvex == NULL)
  204.                ermes("Memory: polygons");
  205.          }
  206.          POLY[npoly++] = i;
  207.       }
  208.       if(npoly == 1)
  209.          ermes("Internal error: npoly = 1");
  210.       if(npoly == 2) {
  211.          add_linesegment(POLY[0], POLY[1]);
  212.          continue;
  213.       }
  214.       jface = j;
  215.       if(!counter_clock(0, 1, 2, 0))
  216.          continue;         /* backface */
  217.       for(i=1; i<=npoly; i++) {
  218.          ii = i%npoly;
  219.          code = POLY[ii];
  220.          vertexnr = abs(code);
  221.          if(code<0)
  222.             POLY[ii] = vertexnr;
  223.          else
  224.             add_linesegment(POLY[i-1], vertexnr);
  225.       }
  226.  
  227.       /* Triangulation of a polygon */
  228.  
  229.       count = 1;
  230.       polygonconvex = 0;
  231.       i1 = -1;
  232.       while(npoly > 2) {
  233.          if(!polygonconvex) {
  234.             polygonconvex = 1;
  235.             for(ii=0; ii<npoly; ii++) {
  236.                i0 = (ii == 0 ? npoly-1 : ii-1);
  237.                i2 = (ii == npoly-1 ? 0 : ii+1);
  238.                vertexconvex[ii] = counter_clock(i0, ii, i2, 0);
  239.                if(!vertexconvex[ii])
  240.                   polygonconvex = 0;
  241.             }
  242.          }
  243.          loopcount = npoly;
  244.          do {
  245.             if(++i1 >=npoly)
  246.                i1 = 0;
  247.                 i0 = (i1 == 0 ? npoly-1 : i1 - 1);
  248.                 i2 = (i1 == npoly-1 ? 0 : i1 + 1);
  249.          }
  250.          while(--loopcount && !polygonconvex &&
  251.               (!vertexconvex[i1] || includesvertex(i0, i1, i2)));
  252.               /* store traingle: */
  253.          counter_clock(i0, i1, i2, count++);
  254.          npoly--;
  255.          for(ii=i1; ii<npoly; ii++)
  256.             POLY[ii] = POLY[ii+1];
  257.       }
  258.    }
  259.    /* Add nearest triangles to screen lists:    */
  260.    for(ipix=0; ipix<NSCREEN; ipix++)
  261.    for(jpix=0; jpix<NSCREEN; jpix++) {
  262.       pointer = &(SCREEN[ipix][jpix]);
  263.       if(pointer->tr_cov >= 0) {
  264.          pnode = TRLIST + iTRLIST;
  265.          pnode->ptria = pointer->tr_cov;
  266.          pnode->next  = pointer->start;
  267.          pointer->start = iTRLIST;
  268.          if(++iTRLIST == nTRLIST)
  269.             ermes("Memory: Triangle list");
  270.       }
  271.    }
  272.    farfree(POLY);
  273.    farfree(vertexconvex);
  274.    ntrsetalloc = 1000;
  275.    trset = (int *)farcalloc(ntrsetalloc, sizeof(unsigned int));
  276.    if(trset == NULL)
  277.       ermes("Memory: set of traingles");
  278.    lsegstart = NULL;
  279.    clearpage();
  280.    init_viewport();
  281.  
  282.    /* Draw all line segments as far as they are visible */
  283.    for(P=1; P<=maxvertex; P++) {
  284.       pvertex = VERTEX + P;   /* = &VERTEX[P]   */
  285.       ptr = pvertex->connect;
  286.       if(ptr == dummy || ptr == NULL)
  287.          continue;
  288.       xP = pvertex->xe;
  289.       yP = pvertex->ye;
  290.       zP = pvertex->ze;
  291.       XP = xP/zP;
  292.       YP = yP/zP;
  293.       for( ; ptr != NULL; ptr = ptr->next) {
  294.             Q = ptr->i;
  295.          pvertex = VERTEX + Q;   /* = &VERTEX[Q]   */
  296.          xQ = pvertex->xe;
  297.          yQ = pvertex->ye;
  298.          zQ = pvertex->ze;
  299.          XQ = xQ/zQ;
  300.          YQ = yQ/zQ;
  301.          /* Using the screen lists, we shall build the
  302.             set of triangles that may hide points of PQ: */
  303.          if(XP<XQ || (XP == XQ && YP < YQ)) {
  304.             Xlft  = XP;
  305.                 Ylft  = YP;
  306.             Xrght = XQ;
  307.             Yrght = YQ;
  308.          }
  309.          else {
  310.             Xlft  = XQ;
  311.                 Ylft  = YQ;
  312.             Xrght = XP;
  313.             Yrght = YP;
  314.          }
  315.          ipixleft  = xwhole(Xlft);
  316.          ipixright = xwhole(Xrght);
  317.             denom = Xrght - Xlft;
  318.          if(ipixleft != ipixright)
  319.             slope = (Yrght - Ylft)/denom;
  320.          jbot = jtop = ywhole(Ylft);
  321.          for(ipix=ipixleft; ipix<=ipixright; ipix++) {
  322.             if(ipix == ipixright)
  323.                jI = ywhole(Yrght);
  324.             else
  325.                jI = ywhole(Ylft+(xreal(ipix-1)-Xlft)*slope);
  326.             LOWER[ipix] = min(jbot, jI);
  327.             jbot = jI;
  328.             UPPER[ipix] = max(jtop, jI);
  329.             jtop = jI;
  330.          }
  331.          ntrset = 0;
  332.          for(ipix=ipixleft; ipix<=ipixright; ipix++)
  333.          for(jpix=LOWER[ipix]; jpix<=UPPER[ipix]; jpix++) {
  334.             pointer = &(SCREEN[ipix][jpix]);
  335.             inode = pointer->start;
  336.             while(inode >= 0) {
  337.             /* At the end of the list we have inode = -1 */
  338.                pnode = TRLIST + inode;
  339.                itria = pnode->ptria;
  340.                /* The triangle will be stored only if it is not
  341.                   yet present in the array trset (the traingle set). */
  342.                if(ntrset == ntrsetalloc) {
  343.                   ntrsetalloc += (ntrsetalloc >> 2);
  344.                   trset = (int *)farrealloc(trset,
  345.                      (long)ntrsetalloc*sizeof(unsigned int));
  346.                   if(trset == NULL)
  347.                      ermes("Memory: traingle set");
  348.                }
  349.                trset[ntrset] = itria;  /* sentinel */
  350.                jtr = 0;
  351.                while(trset[jtr] != itria)
  352.                   jtr++;
  353.                if(jtr == ntrset) {
  354.                   ntrset++;      /* store triangle */
  355.                }
  356.                inode = pnode->next;
  357.             }
  358.          }
  359.          /* Now trset[0], ..., trset[ntrset-1] is the set of
  360.             traingles that may hide points of PQ.           */
  361.  
  362.          linesegment(xP, yP, zP, xQ, yQ, zQ, 0);
  363.          while(lsegstart != NULL) {
  364.             xxP = lsegstart->xP;
  365.             yyP = lsegstart->yP;
  366.             zzP = lsegstart->zP;
  367.             xxQ = lsegstart->xQ;
  368.             yyQ = lsegstart->yQ;
  369.             zzQ = lsegstart->zQ;
  370.             kk  = lsegstart->k;
  371.             auxlseg = lsegstart;
  372.             lsegstart = lsegstart->next;
  373.             farfree(auxlseg);
  374.             linesegment(xxP, yyP, zzP, xxQ, yyQ, zzQ, kk);
  375.          }
  376.       }
  377.    }
  378.    farfree(dummy);
  379.    farfree(trset);
  380.    for(i=1; i<=maxvertex; i++)
  381.    if(VERTEX[i].inuse) {
  382.       dummy = VERTEX[i].connect;
  383.       while(dummy != NULL) {
  384.          dummy0 = dummy;
  385.          dummy = dummy->next;
  386.          farfree(dummy0);
  387.       }
  388.    }
  389.    farfree(TRLIST);
  390.    farfree(TRIANGLE);
  391. }
  392.  
  393. void coeff(float rho, float theta, float phi)
  394. {
  395.    double th, ph, costh, sinth, cosph, sinph;
  396.  
  397.    /* Angles in radians: */
  398.    th = theta * PIdiv180;
  399.    ph = phi   * PIdiv180;
  400.  
  401.    costh = cos(th);
  402.    sinth = sin(th);
  403.    cosph = cos(ph);
  404.    sinph = sin(ph);
  405.  
  406.             /* Elements of viewing matrix V: */
  407.  
  408.    v11 = -sinth;        v12 = -cosph*costh;  v13 = -sinph*costh;
  409.    v21 =  costh;        v22 = -cosph*sinth;  v23 = -sinph*sinth;
  410.    v31 =    0  ;        v32 =  sinph;        v33 = -cosph;
  411.                                              v43 =  rho;
  412.  
  413. }
  414.  
  415. static void add_linesegment(int P, int Q)
  416.  
  417. /* This function may insert line segment PQ in the linked
  418.    list starting in VERTEX[p].connect. (If necessary, we
  419.    first exchange the values of P and Q such that P becomes
  420.    the smaller of the two.) If PQ has been stored in the
  421.    list previously, it is not duplicated. If a previously
  422.    stored segment PQ is found, we compare the face to which
  423.    it belongs with the face to which the new one belongs.
  424.    If the normal vectors to these two faces are almost
  425.    parallel (the inner product between their normal vectors
  426.    being greater than 'surflimit'), the old line segment is
  427.    deleted.
  428. */
  429. {
  430.     struct linsegface *ptr, *new, **pp, **pp0;
  431.    int iaux;
  432.  
  433.    if(P > Q) {
  434.       iaux = P;
  435.       P    = Q;
  436.       Q    = iaux;
  437.    }
  438.    /* Now: P < Q */
  439.  
  440.     pp0 = pp = &(VERTEX[P].connect);
  441.    ptr = *pp;
  442.    while(ptr != NULL && ptr->i != Q) {
  443.       pp = &(ptr->next);
  444.       ptr = *pp;
  445.    }
  446.    if(ptr == NULL) {
  447.       new = lsegfacenode();
  448.       new->i = Q;
  449.       new->a = a;
  450.       new->b = b;
  451.       new->c = c;
  452.       new->next = *pp0;
  453.       *pp0 = new;
  454.    }
  455.    else if (a*(ptr->a) + b*(ptr->b) + c*(ptr->c) > surfl) {
  456.       *pp = ptr->next;
  457.       farfree(ptr);
  458.    }
  459. }
  460.  
  461. static int counter_clock(int i0, int i1, int i2, int code)
  462. /* code = 0: compute orientation;
  463.    code = 1: compute a, b, c, h; store the first triangle;
  464.    code > 1: check if next triangle is coplanar; store it
  465. */
  466. {
  467.     int A=POLY[i0], B=POLY[i1], C=POLY[i2], i, l;
  468.    double xA, yA, zA, xB, yB, zB, xC, yC, zC, r,
  469.           XA, YA, XB, YB, XC, YC, h0,
  470.           DA, DB, DC, D, DAB, DAC, DBC, aux, dist,
  471.           xR, yR, dev, tot;
  472.    static double h;
  473.  
  474.    if(A < 0) A = -A;
  475.    if(B < 0) B = -B;
  476.    if(A < 0) C = -C;
  477.  
  478.    pvertex = VERTEX + A;
  479.    xA = pvertex->xe;
  480.    yA = pvertex->ye;
  481.    zA = pvertex->ze;
  482.  
  483.    pvertex = VERTEX + B;
  484.    xB = pvertex->xe;
  485.    yB = pvertex->ye;
  486.    zB = pvertex->ze;
  487.  
  488.    pvertex = VERTEX + C;
  489.    xC = pvertex->xe;
  490.    yC = pvertex->ye;
  491.    zC = pvertex->ze;
  492.  
  493.    h0 = xA * (yB*zC - yC*zB) -
  494.         xB * (yA*zC - yC*zA) +
  495.         xC * (yA*zB - yB*zA);
  496.  
  497.    if(code == 0) {
  498.       if(h0 <= eps)
  499.          return 0;
  500.       a =   yA*(zB-zC) - yB*(zA-zC) + yC*(zA-zB) ;
  501.       b = -(xA*(zB-zC) - xB*(zA-zC) + xC*(zA-zB));
  502.       c =   xA*(yB-yC) - xB*(yA-yC) + xC*(yA-yB);
  503.       r = sqrt(a*a + b*b + c*c);
  504.       if(r == 0.0)
  505.          r = eps;
  506.       a /= r;
  507.       b /= r;
  508.       c /= r;
  509.       h = h0/r;
  510.       return 1;
  511.    }
  512.  
  513.    /* If h0=0, plane ABC passes through E and hides nothing.
  514.       If h0<0, triangle ABC is a backface.
  515.       In both cases, iTRIANGLE is not incremented and the
  516.       traingles of the polygon are not stored.
  517.       a, b, c will be used in add_linesegment.
  518.    */
  519.  
  520.    dev = fabs(a*xC + b*yC + c*zC - h);
  521.    tot = eps + 0.001*fabs(h);
  522.    if(dev > tot) {
  523.       to_text();
  524.       printf("Vertices not in the same plane:\n");
  525.       for(i=1; i<pface[jface][0]; i++)
  526.          printf("%d ", pface[jface][i]);
  527.       exit(1);
  528.    }
  529.  
  530.    ptriangle = TRIANGLE + iTRIANGLE;
  531.    ptriangle->A = A;
  532.    ptriangle->B = B;
  533.    ptriangle->C = C;
  534.    ptriangle->a = a;
  535.    ptriangle->b = b;
  536.    ptriangle->c = c;
  537.    ptriangle->h = h;
  538.  
  539.    /* The triangle will be now be stored in the screen lists
  540.       of the associated pixels; first the array LOWER,
  541.       UPPER, LOW, UP are defined:
  542.    */
  543.  
  544.    XA = xA/zA;
  545.    YA = yA/zA;
  546.    XB = xB/zB;
  547.    YB = yB/zB;
  548.    XC = xC/zC;
  549.    YC = yC/zC;
  550.  
  551.    DA = XB*YC - XC*YB;
  552.    DB = XC*YA - XA*YC;
  553.    DC = XA*YB - XB*YA;
  554.  
  555.    D  = DA + DB + DC;
  556.  
  557.    DAB= DC - M*(XA-XB);
  558.    DAC= DB - M*(XC-XA);
  559.    DBC= DA - M*(XB-XC);
  560.  
  561.    topcode[0] = (D*DAB>0);
  562.    topcode[1] = (D*DAC>0);
  563.    topcode[2] = (D*DBC>0);
  564.  
  565.    Xleft[0] = XA; Yleft[0] = YA; Xright[0] = XB; Yright[0] = YB;
  566.    Xleft[1] = XA; Yleft[1] = YA; Xright[1] = XC; Yright[1] = YC;
  567.    Xleft[2] = XB; Yleft[2] = YB; Xright[2] = XC; Yright[2] = YC;
  568.  
  569.    for(l=0; l<3; l++)      /* l = triangle - side number */
  570.    if(Xleft[l] > Xright[l] ||
  571.      (Xleft[l] == Xright[l] && Yleft[l] > Yright[l])) {
  572.       aux = Xleft[l]; Xleft[l] = Xright[l]; Xright[l] = aux;
  573.       aux = Yleft[l]; Yleft[l] = Yright[l]; Yright[l] = aux;
  574.    }
  575.    ipixmin = xwhole(min3(XA,XB,XC));
  576.    ipixmax = xwhole(max3(XA,XB,XC));
  577.    for(ipix=ipixmin; ipix<=ipixmax; ipix++) {
  578.       LOWER[ipix] = UP[ipix] = 10000;
  579.       UPPER[ipix] = LOW[ipix] = -10000;
  580.    }
  581.  
  582.    for(l=0; l<3; l++) {
  583.       ipixleft  = xwhole(Xleft[l]);
  584.       ipixright = xwhole(Xright[l]);
  585.       denom = Xright[l] - Xleft[l];
  586.       if(ipixleft != ipixright)
  587.          slope = (Yright[l] - Yleft[l])/denom;
  588.       j_old = ywhole(Yleft[l]);
  589.       for(ipix=ipixleft; ipix<=ipixright; ipix++) {
  590.          if(ipix == ipixright)
  591.             jI = ywhole(Yright[l]);
  592.          else
  593.             jI = ywhole(Yleft[l] + (xreal(ipix+1) - Xleft[l])*slope);
  594.          if(topcode[l]) {
  595.             UPPER[ipix] = max3(j_old, jI, UPPER[ipix]);
  596.             UP[ipix] = min3(j_old, jI, UP[ipix]);
  597.          }
  598.          else {
  599.             LOWER[ipix] = min3(j_old, jI, LOWER[ipix]);
  600.             LOW[ipix] = max3(j_old, jI, LOW[ipix]);
  601.          }j_old = jI;
  602.       }
  603.    }
  604.  
  605.    /* For screen column ipix, the triangle is associated
  606.       only with pixels in the rows
  607.       LOWER[ipix], ..., UPPER[ipix].
  608.       The subrange LOW[ipix]+1,...,UP[ipix]-1 of these rows
  609.       denote pixels that lie completely within the triangle.
  610.    */
  611.  
  612.    for(ipix=ipixmin; ipix<=ipixmax; ipix++)
  613.     for(jpix=LOWER[ipix]; jpix<=UPPER[ipix]; jpix++) {
  614.       pointer = &(SCREEN[ipix][jpix]);
  615.       if(jpix>LOW[ipix] && jpix<UP[ipix]) {
  616.          xR = xmin + (ipix + 0.5)*deltaX;
  617.             yR = ymin + (jpix + 0.5)*deltaY;
  618.          denom = a*xR + b*yR + c*d;
  619.          dist = fabs(denom)>eps ? h*sqrt(xR*xR + yR*yR + 1.)/denom : BIG;
  620.          /* The line drawn from the viewpoint E to pixel point (xR, yR, 1)
  621.             intersects plane ABC at a distance dist from E.
  622.          */
  623.          if(dist < pointer->tr_dist) {
  624.             pointer->tr_cov = iTRIANGLE;
  625.             pointer->tr_dist = dist;
  626.          }
  627.       }
  628.       else {      /* Add triangle to screen list: */
  629.          pnode = TRLIST + iTRLIST;
  630.          pnode->ptria = iTRIANGLE;
  631.          pnode->next = pointer->start;
  632.          pointer->start = iTRLIST;
  633.          if(++iTRLIST == nTRLIST)
  634.             ermes("Memory: Triangle list");
  635.  
  636.       }
  637.    }
  638.     if(++iTRIANGLE == nTRIANGLE) 
  639.         ermes("Memory: Triangles");
  640.    return 1;
  641. }
  642.  
  643. static void viewing(double x, double y, double z,
  644.             double *pxe, double *pye, double *pze)
  645. {
  646.    *pxe = v11*x + v21*y;
  647.    *pye = v12*x + v22*y + v32*z;
  648.    *pze = v13*x + v23*y + v33*z + v43;
  649. }
  650.  
  651.  
  652. static void linesegment(double xP, double yP, double zP,
  653.          double xQ, double yQ, double zQ, int k)
  654.    /* Line segment PQ is to be drawn, as far as it is not
  655.       hidden by triangles trset[k] to trset[ntrset-1].  */
  656. {
  657.    int A, B, C, i, Pbeyond, Qbeyond,
  658.       outside, Poutside, Qoutside, eA, eB, eC, sum;
  659.    double a, b, c, h, hP, hQ, r1, r2, r3,
  660.       xA, yA, zA, xB, yB, zB, xC, yC, zC,
  661.       dA, dB, dC, labmin, labmax, lab, mu,
  662.       xmin, ymin, zmin, xmax, ymax, zmax,
  663.       C1, C2, C3, K1, K2, K3, denom1, denom2,
  664.       Cpos, Ppos, Qpos, aux, eps1;
  665.    struct triangle *ptriangle;
  666.  
  667.    while(k<ntrset) {
  668.       itria = trset[k];
  669.       ptriangle = TRIANGLE + itria;
  670.         a = ptriangle->a;
  671.         b = ptriangle->b;
  672.         c = ptriangle->c;
  673.         h = ptriangle->h;
  674.  
  675. /* Test 1: */
  676.  
  677.       hP = a*xP + b*yP + c*zP;
  678.       hQ = a*xQ + b*yQ + c*zQ;
  679.       eps1 = eps + eps*h;
  680.       if(hP-h <= eps1 && hQ-h <= eps1) {
  681.          k++;
  682.          continue;      /* PQ not behind ABC */
  683.       }
  684.  
  685. /* Test 2: */
  686.  
  687.       K1 = yP*zQ - yQ*zP;
  688.       K2 = zP*xQ - zQ*xP;
  689.       K3 = xP*yQ - xQ*yP;
  690.  
  691.       A = ptriangle->A;
  692.       B = ptriangle->B;
  693.       C = ptriangle->C;
  694.  
  695.       pvertex = VERTEX + A;
  696.       xA = pvertex->xe;
  697.       yA = pvertex->ye;
  698.       zA = pvertex->ze;
  699.  
  700.       pvertex = VERTEX + B;
  701.       xB = pvertex->xe;
  702.       yB = pvertex->ye;
  703.       zB = pvertex->ze;
  704.  
  705.       pvertex = VERTEX + C;
  706.       xC = pvertex->xe;
  707.       yC = pvertex->ye;
  708.       zC = pvertex->ze;
  709.  
  710.       dA = K1*xA + K2*yA + K3*zA;
  711.       dB = K1*xB + K2*yB + K3*zB;
  712.       dC = K1*xC + K2*yC + K3*zC;
  713.  
  714.       /* If dA, dB and dC have the same sign, the vertices
  715.          A, B, C lie at the same side of plane EPQ.      */
  716.  
  717.       eA = dA>eps ? 1 : dA < meps ? -1 : 0;
  718.       eB = dB>eps ? 1 : dB < meps ? -1 : 0;
  719.       eC = dC>eps ? 1 : dC < meps ? -1 : 0;
  720.  
  721.       sum = eA + eB + eC;
  722.       if(abs(sum) >= 2) {
  723.          k++;
  724.          continue;
  725.       }
  726.       /* If this test succeeds, the (infinite) line PQ    
  727.          lies outside the pyramid EABC (or the line and the 
  728.          pyramid have at most one point in common.
  729.          If the test fails, there is a point if intersection */
  730.  
  731. /* Test 3: */
  732.       Poutside = Qoutside = 0;
  733.       labmin = 1.;
  734.       labmax = 0.;
  735.       for(i=0; i<3; i++) {
  736.          C1 = yA*zB - yB*zA;
  737.          C2 = zA*xB - zB*xA;
  738.          C3 = xA*yB - xB*yA;
  739.          /* C1 x + C2 y + C3 z = 0 is plane EAB */
  740.          Cpos = C1*xC + C2*yC + C3*zC;
  741.          Ppos = C1*xP + C2*yP + C3*zP;
  742.          Qpos = C1*xQ + C2*yQ + C3*zQ;
  743.          denom1 = Qpos - Ppos;
  744.          if(Cpos > eps) {
  745.             Pbeyond = Ppos<meps;
  746.             Qbeyond = Qpos<meps;
  747.             outside = Pbeyond && Qpos<=eps || Qbeyond && Ppos<=eps;
  748.          }
  749.          else if(Cpos<meps) {
  750.             Pbeyond = Ppos>eps;
  751.             Qbeyond = Qpos>eps;
  752.             outside = Pbeyond && Qpos>=meps || Qbeyond && Ppos>=meps;
  753.          }
  754.          else 
  755.             outside = 1;
  756.          if(outside)
  757.             break;
  758.          lab = fabs(denom1) <= eps ? 1.e7 : -Ppos/denom1;
  759.          /* lab indicates where PQ meet plane EAB */
  760.          Poutside |= Pbeyond;
  761.          Qoutside |= Qbeyond;
  762.          denom2 = dB - dA;
  763.          mu = fabs(denom2) <= eps ? 1.e7 : -dA/denom2;
  764.          /* mu tells where AB meets plane EPQ */
  765.          if(mu >= meps && mu <= oneplus && lab >= meps && lab <= oneplus) {
  766.             if(lab < labmin)
  767.                labmin = lab;
  768.             if(lab > labmax)
  769.                labmax = lab;
  770.          }
  771.          aux = xA; xA = xB; xB = xC; xC = aux;
  772.          aux = yA; yA = yB; yB = yC; yC = aux;
  773.          aux = zA; zA = zB; zB = zC; zC = aux;
  774.          aux = dA; dA = dB; dB = dC; dC = aux;
  775.       }     
  776.       if(outside) {
  777.          k++;
  778.             continue;
  779.       }
  780.  
  781. /* Test 4: */
  782.       if(!(Poutside || Qoutside))
  783.          return;                 /* PQ invisible */
  784.  
  785. /* Test 5: */
  786.       r1 = xQ - xP;
  787.       r2 = yQ - yP;
  788.       r3 = zQ - zP;
  789.  
  790.       xmin = xP + labmin*r1;
  791.       ymin = yP + labmin*r2;
  792.       zmin = zP + labmin*r3;
  793.       if(a*xmin + b*ymin + c*zmin -h < -eps1) {
  794.          k++;
  795.          continue;
  796.       }
  797.       xmax = xP + labmax*r1;
  798.       ymax = yP + labmax*r2;
  799.       zmax = zP + labmax*r3;
  800.       if(a*xmax + b*ymax + c*zmax -h < -eps1) {
  801.          k++;
  802.          continue;
  803.       }
  804.       /* If this test succeeds, an intersection of PQ
  805.          and the pyramid lies in front of plane ABC.  */
  806.  
  807. /* Test 6: */
  808.       if(Poutside || hP < h - eps1) {
  809.          auxlseg = lsegstart;
  810.             lsegstart = lsegnode();
  811.          lsegstart->xP = xP;
  812.          lsegstart->yP = yP;
  813.          lsegstart->zP = zP;
  814.          lsegstart->xQ = xmin;
  815.          lsegstart->yQ = ymin;
  816.          lsegstart->zQ = zmin;
  817.          lsegstart->k  = k+1;
  818.          lsegstart->next = auxlseg;
  819.       }
  820.       if(Qoutside || hQ < h - eps1) {
  821.          auxlseg = lsegstart;
  822.             lsegstart = lsegnode();
  823.          lsegstart->xP = xmax;
  824.          lsegstart->yP = ymax;
  825.          lsegstart->zP = zmax;
  826.          lsegstart->xQ = xQ;
  827.          lsegstart->yQ = yQ;
  828.          lsegstart->zQ = zQ;
  829.          lsegstart->k  = k+1;
  830.          lsegstart->next = auxlseg;
  831.       }
  832.       return;
  833.    }
  834.  
  835.    move(d*xP/zP+c1, d*yP/zP+c2);
  836.    draw(d*xQ/zQ+c1, d*yQ/zQ+c2);
  837. }
  838.  
  839. void init_viewport(void)
  840. {
  841.    float len = 0.3, len1 = 0.2;
  842.    move(0.0, len);
  843.    draw(0.0, 0.0);
  844.    draw(len, 0.0);
  845.    move(x_max-len, 0.0);
  846.    draw(x_max, 0.0);
  847.    draw(x_max, len);
  848.    move(x_max, y_max-len1);
  849.     draw(x_max, y_max);
  850.    draw(x_max-len1, y_max);
  851.    move(len1, y_max);
  852.    draw(0.0, y_max);
  853.    draw(0.0, y_max-len1);
  854.    /* The distinction between len and len1 enables us to
  855.       tell the top from the bottom                       */
  856.  
  857. }
  858.  
  859. static int includesvertex(int h, int i, int j)
  860. /* Does triangle hij include another vertex? */
  861. {
  862.     int l;
  863.    double zh, zi, zj, zl, XH, YH, XI, YI, XJ, YJ, XL, YL;
  864.  
  865.    pvertex = VERTEX + POLY[h]; zh = pvertex->ze;
  866.    XH = pvertex->xe / zh;  YH = pvertex->ye /zh;
  867.  
  868.    pvertex = VERTEX + POLY[i]; zi = pvertex->ze;
  869.    XI = pvertex->xe / zi;  YI = pvertex->ye /zi;
  870.  
  871.    pvertex = VERTEX + POLY[j]; zj = pvertex->ze;
  872.    XJ = pvertex->xe / zj;  YJ = pvertex->ye /zj;
  873.  
  874.    for (l=0; l<npoly; l++) {
  875.       if(!vertexconvex[l] && l != h && l != i && l != j) {
  876.          pvertex = VERTEX + POLY[l]; zl = pvertex->ze;
  877.          XL = pvertex->xe / zl;
  878.          YL = pvertex->ye / zl;
  879.  
  880.          if(sameside(XH, YH, XL, YL, XI, YI, XJ, YJ) &&
  881.             sameside(XI, YI, XL, YL, XJ, YJ, XH, YH) &&
  882.             sameside(XJ, YJ, XL, YL, XH, YH, XI, YI))
  883.           return 1;
  884.       }
  885.    }
  886.    return 0;
  887. }
  888.  
  889. static int sameside(double xj, double yj, double xl, double yl,
  890.                     double xh, double yh, double xi, double yi)
  891.  
  892. /* Do points j and l lie on the same side of line hi? */
  893. {
  894.    double a, b, c;
  895.  
  896.    a = yh - yi;
  897.    b = xi - xh;
  898.    c = xh*yi - yh*xi;
  899.    return (a*xj + b*yj + c) * (a*xl + b*yl + c) > 0.0;
  900. }
  901.  
  902. struct lseglist *lsegnode(void)
  903. {
  904.    struct lseglist *p;
  905.  
  906.    p = (struct lseglist *)
  907.        farmalloc((long)sizeof(struct lseglist));
  908.    if(p == NULL)
  909.       ermes("Memory: line segment PQ");
  910.    return p;
  911. }
  912.  
  913. struct linsegface *lsegfacenode(void)
  914. {
  915.    struct linsegface *p;
  916.  
  917.    p = (struct linsegface *)
  918.        farmalloc((long)sizeof(struct linsegface));
  919.    if(p == NULL)
  920.       ermes("Memory: edges");
  921.    return p;
  922. }
  923.  
  924. static int allocsize(int divfactor, int elsize)
  925. {
  926.    long lll;
  927.  
  928.    lll = (farcoreleft()/divfactor)/elsize;
  929.    return (lll < 32767 ? (int)lll : 32767);
  930. }
  931.  
  932.  
  933.  
  934.  
  935.  
  936.  
  937.  
  938.